/***************************************************************
                       ends2ends
Routines Ends->Ends
**************************************************************/
#include <stdio.h>
#include <malloc.h>
#include <pwd.h>
#include <limits.h>

#ifndef MAC_BUILD
#include <wait.h>
#else
#include <sys/wait.h>
#endif

#include "clam.h"
#include <gtk/gtkwidget.h>
#include <sys/times.h>
#include "seq.h"

#define Show {if (graphActivate(g92)){graphPop(); graphRedraw();}}
#define Output {if (Pz.logFlag) {fprintf(Zlogfp,"%s",ZBuf); fflush(Zlogfp);} \
      if (!Zbatch_flag && Pz.stdFlag) {printf("%s",ZBuf); fflush(stdout);}}

extern int g92;
extern int Num_threads;
extern int ZincrFlag;
extern int called_by_merge;
extern char scoreText[60];
extern int twoEndsFlag;
extern FILE* graphQueryOpen();
extern int loadCz(),  ZscoreCpM(), ZboolCpM();
extern void move_selected(), fix_cb_offsets(), Zproj_results();
extern void scanClam2(), clearMsg(), AutoCtgMsg();
extern void precompute_sulston();
extern void pr_end_timer(char *msg, clock_t real, struct tms * tmsstart, struct tms *tmssend);
extern void Zfree_cb_data_structures();
static double find_min_eval(double* scorelist1,double* scorelist2, int count1, int count2);
static double find_nth_eval(double* scorelist,int count, int n);
static int eval_cmp(const void *pd1, const void *pd2);
extern void redo_sequence_placements();
extern int interval_gap(int s1,int e1,int s2,int e2);

void append_remark(int clone, char* rem);

static void update_end_to_end_struc ( int *);
static void update_ZBufPtr();
void sig_hdlr();
void flipcontig(int ctg);
static void flip_tctgs(int ctg);

char * ZBufPtr = NULL;
unsigned int ZBufSize = 4096;
char * ZBufPrint = NULL;
int ends_ends_seq_confirm = 0;
extern GtkWidget *ctg_window;

#define END 75          /* max number of end clones we'll use from a contig */
#define MERGE_MAX 10    /* used to make the match array large enough (hopefully) */

static struct tmpctg
{
   int n;           /* number of end clones this contig has */
   int cin[END];    /* indices of the end clones */
   char w[END];     /* types of the end clones */
   int merge_to;    /* contig merged to, for following merge chains */
   int num_merges;  /* total number of merges for this contig */
   int nL;          /* number of L clones */
   int nR;          /* number of R clones */
   int nB;          /* number of B clones */
   int redoCB;      /* needed for the file-driven merging */
   int Rmerge;      /* ctg merged to the R side */
   int Lmerge;      /* ctg merged to the L side */
   int flipped;     /* ctg has been flipped in a merge */
} *tctg;

static struct match_info 
{
   int ctg1;
   int ctg2;
   char type1;
   char type2;
   int matches;
   int clones1;
   int clones2;
   int ctgnum1;
   int ctgnum2;
   int best1;
   int best2;
   double best_score;
   double min_eval;
   int cancelled;
} *match_list, *next_match_entry, *last_match_entry;

/* the thread structure for end-end parallelizing */ 
struct end_to_end_threads
{
  int start_contig;
  int end_contig;
}end_to_end[500];

int matches1[END]; /* these store left/right matched end clones*/
int matches2[END];

FPC_TREE ctg_pair_tree;
int seq_pair_tree_filled = 0;

// this is sort of a hack but it provides the needed speed
int merge_confirmed_by_sequence(int ctg1, int ctg2)
{
    FPC_SEQ* pseq;
    FPC_TREE* pnode;
    FPC_STACK tempstack;
    struct seqctgpos *ploc1, *ploc2;
    char temp[100];
    int e1, e2, cspan, sspan;

    /* build the tree of pairs indicated by sequence */
    if (!seq_pair_tree_filled)
    {
        fpc_stack_init(&tempstack);  
        pnode = &seq_tree;

        while ((pnode = fpc_tree_next_node(pnode, &tempstack))) 
        {     
            if (!pnode->pdata) continue;
            pseq = (FPC_SEQ*)pnode->pdata;
            if (pseq->parent) continue;

            for (ploc1 = pseq->ctgpos; ploc1; ploc1 = ploc1->next)
            {
                ctg1 = ploc1->ctg;
                if (contigs[ctg1].ctgstat == DEAD || contigs[ctg1].ctgstat == AVOID)
                {
                    continue;
                }
                e1 = ploc1->ctg_start;
                e2 = contigs[ctg1].right - ploc1->ctg_end;
                cspan = ploc1->ctg_end - ploc1->ctg_start;
                sspan = ploc1->seq_end - ploc1->seq_start;
                if (0) //cspan < Pz.minspanInt || sspan < Pz.seq_minspanInt )
                {
                    continue;
                }

                if (e1 <  Pz.fromendInt  ||  e2 < Pz.fromendInt)
                {
                    // this contig has an endpoint on the sequence. Look for a paired contig. 

                    for (ploc2 = ploc1->next; ploc2; ploc2 = ploc2->next)
                    {
                        ctg2 = ploc2->ctg;
                        if (contigs[ctg2].ctgstat == DEAD || contigs[ctg2].ctgstat == AVOID)
                        {
                            continue;
                        }
                        e1 = ploc2->ctg_start;
                        e2 = contigs[ctg2].right - ploc2->ctg_end;
                        cspan = ploc2->ctg_end - ploc2->ctg_start;
                        sspan = ploc2->seq_end - ploc2->seq_start;
                        if (0) //cspan < Pz.minspanInt || sspan < Pz.seq_minspanInt )
                        {
                            continue;
                        }

                        if (e1 <  Pz.fromendInt  ||  e2 < Pz.fromendInt)
                        {
                            // this contig has an endpoint on the sequence.
                            // see if it's close enough to the previous.

                            if (interval_gap(ploc1->seq_start,ploc1->seq_end,ploc2->seq_start,ploc2->seq_end) <= Pz.seq_fromendInt)
                            {
                                if (ploc1->ctg < ploc2->ctg) sprintf(temp,"%d_%d",ploc1->ctg,ploc2->ctg);
                                else sprintf(temp,"%d_%d",ploc2->ctg,ploc1->ctg);

                                fpc_tree_insert_node(&ctg_pair_tree,(void*)1,temp);

                            }
                        } 
                    }


                } 
            }
		}

        fpc_stack_destroy(&tempstack);
        seq_pair_tree_filled = 1;
    }

        
    if (ctg1 <= ctg2) sprintf(temp,"%d_%d",ctg1,ctg2);
    else sprintf(temp,"%d_%d",ctg2,ctg1);

    if (fpc_tree_find_data(&ctg_pair_tree,temp)) return 1;
    return 0;
}

/***********************************************************
                   DEF: Out_ends
***********************************************************/
void Out_ends(char o1, char o2) /* also used by clam.c for --> Ends */
{
int Zcnt_shared_markers();
    int  mkcnt;
    char  mk[40];

     mkcnt = Zcnt_shared_markers(&C1z,&C2z,&Sz);
     if (mkcnt>0) sprintf(mk,"%2d",mkcnt);
     else mk[0]='\0';

     sprintf(ZBufPtr,               /* CLONE_SZ */
      "Ctg%-4d %c %15s %3db  Ctg%-4d %c %15s %3db  Match %3d %.0e %s\n",
       C1z.ctg, o1, C1z.clone, C1z.nbands,
       C2z.ctg, o2, C2z.clone, C2z.nbands, Sz.match, Sz.prob, mk); 
     update_ZBufPtr();
}
/**************************************************************/
static void update_ZBufPtr()
{
   // Since we don't write more than 100 characters in one go...
   // we subtract - 100 from ZBufSize 
   // This is clumsy but no other option

   if (strlen (ZBufPrint) >= (ZBufSize - 100) )
   {
      ZBufPrint = (char *)realloc (ZBufPrint, ZBufSize * 2);
      ZBufSize *= 2;
   }   
   ZBufPtr = & ZBufPrint[strlen(ZBufPrint)] ;
}

/****************************************************/
static void merge_contigs_by_side(int to, int from, char to_side, char from_side)
{
   int c, i, Q;
   CLONE* clp;
   int offset = 0;

   if (to_side == 'B' || from_side == 'B')
   {
      fprintf(stderr, "FPC Error: B type merge %d%c %d%c\n",to,to_side,from,from_side);
      return;
   }

   /* flip one contig if needed */
   merging = 1;
  if (to_side == from_side)
  {
      if (from_side == 'R')
      {
         sprintf(ZBuf,"Flip contig %d\n",from);Output;             
         flipcontig(from);
         flip_tctgs(from);
         from_side = 'L';
      }
      else
      {
         sprintf(ZBuf,"Flip contig %d\n",to);Output;             
         flipcontig(to);
         flip_tctgs(to);
         to_side = 'R';
      }
  }
   merging = 0;
  
   
   if (to_side == 'R')
   {
      offset = -contigs[from].left + contigs[to].right - Pz.fromendInt;
   }
   else if (from_side == 'R')
   {
      offset = contigs[to].left -contigs[from].right + Pz.fromendInt;
   }
   else
   {
      fprintf(stderr, "FPC Error: confusion in contig merge %d to %d\n",from,to);
      return;
   }

   c =  contigs[from].start;
   for (i = 0; i < contigs[from].count && c != -1; i++)
   {
       clp = arrp(acedata, c, CLONE);
       clp->selected = TRUE;
       c= clp->next;
   } 

   ZincrFlag = 1; /* so move_selected will not change oldctg */
   called_by_merge = 1; PRTMESS=0; 
   Q = contigs[to].ctgQs + contigs[from].ctgQs; /* because move_selected changes the Qs */
   move_selected(from, to, offset); 
   contigs[to].ctgQs = Q;
   contigs[to].approxQs = 1;
   called_by_merge = 0; PRTMESS=1; 
   ZincrFlag = 0; 

   c =  contigs[to].start;
   for (i = 0; i < contigs[to].count && c != -1; i++)
   {
       clp = arrp(acedata, c, CLONE);
       clp->selected = FALSE;
       c= clp->next;
   } 
}

/****************************************************************/
static void init_tctg()
{
  int k, c, n;
  CLONE *clp;
  int left, right;

/* initialize the tmp contig struct, in particular finding the L/R/B clones */
  for (k = 1; k <= max_contig; k++) 
  {
       tctg[k].n = 0;
       tctg[k].merge_to = 0;
       tctg[k].num_merges = 0;
       tctg[k].nL = 0;
       tctg[k].nR = 0;
       tctg[k].nB = 0;       
       tctg[k].Rmerge = 0;       
       tctg[k].Lmerge = 0;       
       tctg[k].redoCB = 0;
       tctg[k].flipped = 0;
  
       if (contigs[k].count == 0) continue;
       if (contigs[k].ctgstat== DEAD || contigs[k].ctgstat == AVOID) continue; 

       c= contigs[k].start;
       while (1)
       {
          clp = arrp(acedata, c, CLONE);
          left = clp->x - contigs[k].left;
          right = contigs[k].right - clp->y;

          if (left <= Pz.fromendInt || 
              right <= Pz.fromendInt )
          {
             n = tctg[k].n;
             if (n >= END) 
             {
                printf("Warning: more than %d end clones for ctg%-3d. Ignoring the rest...\n",
                             END, k);
                break;
             }
             tctg[k].n++;
             tctg[k].cin[n] = c;

			// WMN removed the last test here which was preventing valid merges of small contigs. The 
			// original point of it was no longer so clear either. 			
            // if (left <= Pz.fromendInt && right <= Pz.fromendInt && contigs[k].count <= 2*Pz.mergeMatches )
             if (left <= Pz.fromendInt && right <= Pz.fromendInt )
             {
                tctg[k].w[n] = 'B';
                tctg[k].nB++;
                tctg[k].nL++;
                tctg[k].nR++;
             }
             else if (left < right) 
             {
                tctg[k].w[n] = 'L';
                tctg[k].nL++;
             }
             else 
             {
                tctg[k].w[n] = 'R';
                tctg[k].nR++;
             }
          }
          if (c == contigs[k].last) break;
          else c= clp->next;
       }
   }

}
/***************************************************************
                    add_match
 Add the match info to the match list.
 Also with an auxiliary function to sort by contig num, 
 which helps in later converting RL matches to B. 
*************************************************************/
static void add_match(int i1, int i2, char type1, char type2, 
                      int count1, int count2, int total_count,
                       int best1, int best2, double best_score,
                        double min_eval)
{
    if (i1 >= i2)
    {
       fprintf(stderr,"FPC Error: match pair %d %d not ordered\n",i1,i2);
       return;
    }
    if (best1 > 0 && best2 > 0)
    {
        /* otherwise we came from merge_from_list */
        sprintf(ZBufPtr,"Match: %d%c %d%c cutoff:%.0e\n",contigs[i1].ctg,type1,contigs[i2].ctg,type2,min_eval);
        update_ZBufPtr();
    }
    if (next_match_entry < last_match_entry)
    {
       next_match_entry->ctg1 = i1;      
       next_match_entry->ctg2 = i2;      
       next_match_entry->type1 = type1;      
       next_match_entry->type2 = type2;      
       next_match_entry->clones1 = count1;      
       next_match_entry->clones2 = count2;      
       next_match_entry->matches = total_count;   
       next_match_entry->ctgnum1 = contigs[i1].ctg;      
       next_match_entry->ctgnum2 = contigs[i2].ctg;       
       next_match_entry->best1 = best1;       
       next_match_entry->best2 = best2;       
       next_match_entry->best_score = best_score;       
       next_match_entry->min_eval = min_eval;       
       next_match_entry->cancelled = 0;       

       next_match_entry++;
       tctg[i1].num_merges++;
       tctg[i2].num_merges++;
    }
    else
    {
       sprintf(ZBuf,"Maximum number of matches exceeded\n");Output;
    }
}
/****************************************************************/
static void sort_add_match(int i1, int i2, char type1, char type2, 
                           int count1, int count2, int total_count,
                           int best1, int best2, double best_score,
                            double min_eval)
{
   if (i1 < i2)
   {
      add_match(i1, i2, type1, type2, count1, count2, total_count,
                  best1, best2, best_score, min_eval);
   }
   else if (i2 < i1)
   {
      add_match(i2, i1, type2, type1,count2, count1, total_count,
                best2, best1, best_score, min_eval);
   }
}
/*******************************************************************8
                     check_add_clone
 add the number to the list, if it isn't there
*********************************************************************/
static void check_add_clone(int* list, int* pairlist, double* scorelist, int* count, int c, int c2, double score)
{
   int found = 0;
   int i;

   for (i = 0; i < *count; i++)
   {
      if (list[i] == c)
      {
         found = 1;
         if (score < scorelist[i])
         {
            list[i] = c;
            pairlist[i] = c2;
            scorelist[i] = score;
         }
         break;
      }
   }
   if (found == 0)
   {
      list[*count] = c;
      pairlist[*count] = c2;
      scorelist[*count] = score;
      (*count)++;
   }
   
}
/****************************************************8
                  check_ctg_matches
   Look for matches in the various categories RR etc.. 
   It goes category by category since otherwise it would need 
   a pair of auxiliary arrays like matches1, matches2,
   for each category, so it can track the unique matches for each. 
********************************************************/
static void check_ctg_matches(int i1, int i2)
{
   int n1, n2;
   char types1[] = "RL";
   char types2[] = "RL";
   char* type1;
   char* type2;
   int clones1[END]; /* matching clones on 1st ctg */
   int pairs1[END];  /* best matches of the clones in clones1 */
   double scores1[END]; /* scores of the matches in pairs1 */
   int clones2[END];
   int pairs2[END];
   double scores2[END];
   int count1;
   int count2;
   int total_count;
   int score;
   int best1 = 0;
   int best2 = 0;
   double best_score;
   double mineval;


  if (ends_ends_seq_confirm && !merge_confirmed_by_sequence(i1,i2))
  {
     //sprintf(ZBuf,"\Skipping %d to %d because of lack of sequence confirmation\n",i1,i2); 
     //Output; 
     return;          
  }  


   /* since we have clones which can be both L/R, we save their overlaps
     in order not to repeat them */
   double already_seen[END][END];

   for (n1 = 0; n1 < END; n1++)
   {
      for (n2 = 0; n2 < END; n2++)
      {
         already_seen[n1][n2] = -1;
      }
   }


   for (type1 = types1; *type1 != 0; type1++)
   {
      /* no point looking at this category if there aren't enough clones in it */
      if (  (*type1 == 'R' && tctg[i1].nR < Pz.mergeMatches) ||
             (*type1 == 'L' && tctg[i1].nL < Pz.mergeMatches) )
      {
         continue;
      }

      for (type2 = types2; *type2 != 0; type2++)
      {
         if (  (*type2 == 'R' && tctg[i2].nR < Pz.mergeMatches) ||
                (*type2 == 'L' && tctg[i2].nL < Pz.mergeMatches) ) 
         {
            continue;
         }
         count1 = 0;
         count2 = 0;
         total_count = 0;   
         best_score = 1;

         for (n1 = 0; n1 < tctg[i1].n; n1++)
         {
             if (tctg[i1].w[n1] == *type1 || tctg[i1].w[n1] == 'B')
             {
                loadCz(&C1z, tctg[i1].cin[n1]);
                for (n2 = 0; n2 < tctg[i2].n; n2++)
                {
                   if (tctg[i2].w[n2] == *type2 || tctg[i2].w[n2] == 'B')
                   {
                      loadCz(&C2z, tctg[i2].cin[n2]);
                      if (already_seen[n1][n2]  != -1)
                      {
                         Sz.prob = already_seen[n1][n2];
                      }
                      else
                      {
                         Zsulston(&C1z,&C2z,&Sz);
                      }
                      score = ZscoreCpM(&C1z,&C2z,&Sz);               
                      if (score > 0) 
                      {
                         check_add_clone(clones1,pairs1,scores1,&count1,tctg[i1].cin[n1],tctg[i2].cin[n2],Sz.prob);
                         check_add_clone(clones2,pairs2,scores2,&count2,tctg[i2].cin[n2],tctg[i1].cin[n1],Sz.prob);
                         total_count++;
                         if (already_seen[n1][n2] == -1)
                         {
                            Out_ends(tctg[i1].w[n1], tctg[i2].w[n2]);
                         }
                         
                         if (Sz.prob < best_score)
                         {
                            best_score = Sz.prob;
                            best1 = tctg[i1].cin[n1];
                            best2 = tctg[i2].cin[n2];
                         }
                      }
                      already_seen[n1][n2] = Sz.prob;
                   }
                } 
             }
         }
         if (count1 >= Pz.mergeMatches && count2 >= Pz.mergeMatches)
         {
            mineval = find_min_eval(scores1,scores2, count1, count2);
            sort_add_match(i1,i2,*type1,*type2,count1,count2,total_count,
                                best1, best2, best_score,mineval);
         }
      }
   }
}

/********************************************************************
  Find the minimum e-value for which there are enough overlaps for 
    the merge to happen. It's not strictly correct because only the
    best matches for each clone are stored and this could affect the 
    result, but probably not much.
********************************************************************/
static double find_min_eval(double* scorelist1,double* scorelist2, int count1, int count2)
{
    double e1;
    double e2;

    e1 = find_nth_eval(scorelist1, count1, Pz.mergeMatches);
    e2 = find_nth_eval(scorelist2, count2, Pz.mergeMatches);

    return (e1 > e2 ? e1 : e2);

}
static double find_nth_eval(double* scorelist,int count, int n)
{
    double *sorted_scores;
    double ret;

    if (count < n)
    {
        fprintf(stderr,"FPC Error scoring merge:count %d, n %d\n",count,n);
        return 0.0;
    }

    sorted_scores = (double*)malloc(count*sizeof(double));
    memcpy(sorted_scores,scorelist,count*sizeof(double));
    qsort(sorted_scores,count, sizeof(double), eval_cmp);   
    
    ret = sorted_scores[Pz.mergeMatches - 1];
    free(sorted_scores);
    return ret;
}
/* Compare e-values for qsort. Lower evalue will get sorted to the front.  */
static int eval_cmp(const void *pd1, const void *pd2)
{
    double d1, d2;
    d1 = *(double*)pd1;
    d2 = *(double*)pd2;
    if (d1 > d2)
    {
        return 1;
    }
    else if (d1 < d2)
    {
        return -1;
    }
    return 0;
}
/********************************************************************/
static void safe_strcat(char* dest, char* new, int max)
{
   if (strlen(dest) >= max - 1)
   {
      return;
   }
   strncat(dest, new, max - strlen(dest) - 1);
}

/********************************************************************/
void add_to_projmsg(int k, char* msg)
{
   safe_strcat(contigs[k].projmsg, msg, MSGSIZE);
}

/********************************************************************/
/* find out where a contig has been merged to */
static int locate_contig(int i)
{
    int ret = i;
    while (tctg[ret].merge_to != 0)
    {
       ret = tctg[ret].merge_to;
    }
    return ret;
}
/* flip all the contigs which were merged into the given contig */
static void flip_tctgs(int ctg)
{
   int i;
   for (i=1; i <= max_contig; i++) 
   {
      if (locate_contig(i) == ctg)
      {
         tctg[i].flipped = (tctg[i].flipped == 1 ? 0 : 1);
      }
   }   
}
/*******************************************************************
                       add_match_remarks
   Make contig remarks indicating matches that were found but not carried out.
   This is used for manual merge and for auto, for those contigs with too
   many matches found. 
**********************************************************************/
static void add_match_remarks(int only_for_excess)
{
   struct match_info* match;
   int k;
   char   tmp[MSGSIZE+20]; 

   for (k = 1; k <= max_contig; k++) 
   {
       if (only_for_excess == 1)
       {
          if (tctg[k].num_merges <= MERGE_MAX)
          {
             continue;
          }
          else
          {
             add_to_projmsg(k,"Too Many:");
          }
       }

       /* Contig k is now known to have not been merged to anything.
          The second contig in a pair with k may have moved however */
       for (match = match_list; match < next_match_entry; match++)
       {
          if (match->ctg1 == k)
          {
             sprintf(tmp, " %c%c-%d ctg%d",match->type1,match->type2, 
                        match->matches,match->ctg2);  
             add_to_projmsg(k, tmp);
          
          }
          else if (match->ctg2 == k)
          {
             sprintf(tmp, " %c%c-%d ctg%d",match->type2,match->type1, 
                        match->matches,match->ctg1);   
             add_to_projmsg(k, tmp);
          }
       }
   }
}
/****************************************************************
if the phrase End-merge gets changed here, also change in cb_ok_contigs.c
as it is removed when a new assembled is okay.
*************************************************************/
void add_merge_remarks(int best1, int best2, double best_score)
{
    char rem[COMMENT_SZ];
    CLONE* clp;

    clp = arrp(acedata, best2, CLONE);
    sprintf(rem, "End-merge %.0e %s", best_score, clp->clone);
    append_remark(best1, rem);

    clp = arrp(acedata, best1, CLONE);
    sprintf(rem, "End-merge %.0e %s", best_score, clp->clone);
    append_remark(best2, rem);
}
/********************************************************/
void append_remark(int clone, char* rem)
{
    struct remark* pnew;
    struct remark* pnext;
    CLONE* clp;
    clp = arrp(acedata, clone, CLONE);
    
    pnew = (struct remark* )malloc(sizeof(struct remark));
    strncpy(pnew->message, rem, COMMENT_SZ - 1);
    pnew->box = 1;   /* WN - these don't seem to make a difference */
    pnew->colour = 1;
    pnew->next = 0;

    if (clp->fp_remark == 0)
    {
        clp->fp_remark = pnew;
    }
    else
    {
        for(pnext = clp->fp_remark; pnext->next != 0; pnext = pnext->next);
        pnext->next = pnew;
    }
}


/*************************************************************
    Do the merges, in order of their score
**************************************************************/
static int match_cmp(const void *p1, const void *p2)
{
    struct match_info *pm1, *pm2;

    pm1 = (struct match_info*)p1;   
    pm2 = (struct match_info*)p2;   
    return eval_cmp(&(pm1->min_eval), &(pm2->min_eval));
}
/*********************************************************/
static char reverse_type(char type)
{
    if (type == 'B')
    {
        return type;
    }
    return (type == 'R' ? 'L' : 'R');
}
/*********************************************************/
static char get_side(char type, int ctg)
{
    if (tctg[ctg].flipped == 1)
    {
       return reverse_type(type);
    }
    else return type;
}

static void do_automerges()
{
   struct match_info *match;
   int from, to, k, t;
   char to_side, from_side, c;
   char tmp[MSGSIZE+20];
   int cnt = 0;
   void AutoCtgMsg();

   if(ctg_window!=NULL) gtk_widget_destroy(ctg_window);

   qsort(match_list, (int)(next_match_entry - match_list), sizeof(struct match_info), match_cmp);

   for (k = 1; k <= max_contig; k++) 
   {
       contigs[k].projmsg[0] = 0;
   }

   for (match = match_list; match < next_match_entry; match++)
   {
      if (match->cancelled == 1)
      {
         continue;
      }
      /* check if the end we want to use has been used already */
      if (match->type1 == 'R' && tctg[match->ctg1].Rmerge != 0)
      {
         sprintf(ZBuf,"\nIgnoring merge %d%c to %d%c because of previous merge of %d%c to %d\n",  
                  match->ctg1,match->type1,match->ctg2,match->type2,match->ctg1,match->type1,tctg[match->ctg1].Rmerge); 
         Output; 
         continue;          
      }
      if (match->type1 == 'L' && tctg[match->ctg1].Lmerge != 0)
      {
         sprintf(ZBuf,"\nIgnoring merge %d%c to %d%c because of previous merge of %d%c to %d\n",  
                  match->ctg1,match->type1,match->ctg2,match->type2,match->ctg1,match->type1,tctg[match->ctg1].Lmerge); 
         Output; 
         continue;          
      }
      if (match->type2 == 'R' && tctg[match->ctg2].Rmerge != 0)
      {
         sprintf(ZBuf,"\nIgnoring merge %d%c to %d%c because of previous merge of %d%c to %d\n",  
                  match->ctg1,match->type1,match->ctg2,match->type2,match->ctg2,match->type2,tctg[match->ctg2].Rmerge); 
         Output; 
         continue;          
      }
      if (match->type2 == 'L' && tctg[match->ctg2].Lmerge != 0)
      {
         sprintf(ZBuf,"\nIgnoring merge %d%c to %d%c because of previous merge of %d%c to %d\n",  
                  match->ctg1,match->type1,match->ctg2,match->type2,match->ctg2,match->type2,tctg[match->ctg2].Lmerge); 
         Output; 
         continue;          
      }

      /* The merge is allowed - mark the ends as used */

     switch (match->type1)
     {
        case 'R':
           tctg[match->ctg1].Rmerge = match->ctg2; 
           break;  
        case 'L':
           tctg[match->ctg1].Lmerge = match->ctg2; 
           break; 

     } 

     switch (match->type2)
     {
        case 'R':
           tctg[match->ctg2].Rmerge = match->ctg1; 
           break;  
        case 'L':
           tctg[match->ctg2].Lmerge = match->ctg1; 
           break; 

     } 


      /* One or both of these contigs may have already been merged to another one */
      from = locate_contig(match->ctg1);   
      to = locate_contig(match->ctg2); 
      from_side = get_side(match->type1,match->ctg1);
      to_side = get_side(match->type2,match->ctg2);

      /* Always merge to smaller number contig */

      if (from < to)
      {
         t = from;
         from = to;
         to = t;

         c = from_side;
         from_side = to_side;
         to_side = c;
      }  

      if ( from != to )
      {   
         sprintf(tmp,"End-merge Ctg%d. ",from);
         add_to_projmsg(to,tmp);  
         AutoCtgMsg(to,tmp);
         add_to_projmsg(to,contigs[from].projmsg);  

         sprintf(ZBuf,"\nMerge contig %d%c to %d%c (Original:%d%c %d%c) Score:%.0e\n",  
             from, from_side, to,to_side,match->ctg1,match->type1,match->ctg2,match->type2,match->min_eval); Output; 
         merge_contigs_by_side(to,from,to_side,from_side);

         tctg[from].merge_to = to;  
         tctg[to].redoCB = 1;
 

         add_merge_remarks(match->best1, match->best2, match->best_score);
         cnt++;
      }
      else if (from == to)
      {
          sprintf(ZBuf,"\nIgnoring merge of contig %d to %d (Original:%d%c %d%c)\n",  
             contigs[from].ctg, contigs[to].ctg,match->ctg1,match->type1,match->ctg2,match->type2); 
          Output;
      }

   }


   sprintf(ZBuf,"\nComplete merge of %d contigs.",cnt);Output;
   sprintf(scoreText,"Ends -> Ends  %d merges", cnt); Show;
   Zfree_cb_data_structures();

    redo_sequence_placements();
}



/******************* Parallel Code ***************************************/
/********************************************************************************
       update_end_to_end_struc:

   This functions assigns number of comparisions to be performed per thread 
********************************************************************************/
static void update_end_to_end_struc(int *num_thread)
{
   unsigned long sum=0, num_tries = 0;
   float tries_iterator=0;
   int prev_thread=0;
   int prev_thread1=0;
   float iter_per_thread;
   int i,j,index;
   int actual_num_thread;
   struct end_thread
   {
      int start_contig;
      int end_contig;
   }end_end_thread[500];
   float *cum_tries_contig;

   cum_tries_contig = malloc(sizeof(int)*(max_contig+1));

   if(*num_thread>max_contig)
        *num_thread=max_contig;

   for(i=0;i<*num_thread;i++)
   {
      end_end_thread[i].start_contig=
          end_end_thread[i].end_contig=-1;
   }
   for(i= 0;i<=max_contig; i++)
   {
        cum_tries_contig[i]=0;
   }

   //A cumulative arrray is maintained for the number of comparisions */
   //fill up the cumulative array
   for(i=1 ; i<=max_contig ; i++)
   {
       if( tctg[i].n == 0)
       {
          cum_tries_contig[i] = cum_tries_contig [i-1];
          continue;
       }
       sum=0;
       for(j=i+1;j<=max_contig ; j++)
       {
         sum += tctg[j].n;
       }
       num_tries = tctg[i].n*sum;
       cum_tries_contig[i] = cum_tries_contig[i-1] + num_tries;
   }

       // cum_tries_contig [max_contig] will hold the total num of tries
   iter_per_thread = cum_tries_contig[max_contig] / (*num_thread);
   tries_iterator= iter_per_thread;
   j = 1;
   end_end_thread[0].start_contig= 1;

   //Calculate start and end contig for every thread
   for(index=0; index < *num_thread ; index ++)
   {
        for(i=j ; i<= max_contig; i++)
        {
            if (tctg[i].n == 0)
               continue;
            if(cum_tries_contig[i] >= tries_iterator)
               break;
        }
        if (tries_iterator >= cum_tries_contig [max_contig]) 
        {
              // This is the last thread 
             end_end_thread[index].end_contig = max_contig;
             break;
        }

        if( cum_tries_contig[i] == tries_iterator )
        {
            end_end_thread[index].end_contig = i;
            end_end_thread[index+1].start_contig = i+1;
            j = i+1;
        }
             // Check if > midpt. Also check what the prev thread has selected 
             // if the contig number selected
             // for this iteration is same as the one assigned by the prev thread
        else if ( (tries_iterator - cum_tries_contig [i-1]  > 
                   cum_tries_contig [i] - tries_iterator) || 
                  ((prev_thread || prev_thread1) && (i==j)))
        {
                // Include this contig also
                prev_thread=0;
                prev_thread1=1;
                end_end_thread[index].end_contig = i;
                end_end_thread[index+1].start_contig = i+1;
                j = i+1;
        }
        else
        {
                // Don't include current contig
                prev_thread=1;
                prev_thread1=0;
                end_end_thread[index].end_contig = (i-1) ;
                end_end_thread[index+1].start_contig = i;
                j=i;
        }
        tries_iterator += iter_per_thread;
  }

  /* the max_contig will not participate in any end-end comparisions */
  end_end_thread[index].end_contig = max_contig-1;

  actual_num_thread=-1;

  for(i=0;i<*num_thread;i++)
  {
     if(end_end_thread[i].start_contig==-1 ||
        end_end_thread[i].end_contig == -1 ||
        end_end_thread[i].start_contig > end_end_thread[i].end_contig)
     {
         continue;
     }
     actual_num_thread++;      
     memcpy(&end_to_end[actual_num_thread],&end_end_thread[i],sizeof(struct end_thread));
  }
  *num_thread = actual_num_thread + 1;

  printf("\n Num of threads=%d\n",*num_thread);
  for(i=0;i<*num_thread;i++)
  {
     float num_of_iter = cum_tries_contig [end_to_end[i].end_contig] - 
                                    cum_tries_contig[end_to_end[i].start_contig];
     printf("Thread number %d, start contig %5d, end_contig %5d, Num of computations %.0f \n",
         i,end_to_end[i].start_contig,end_to_end[i].end_contig, num_of_iter);
      fflush(0);
  }
  printf("\n");
  fflush(0);
  free(cum_tries_contig);
}


/****************************************************************
   sig_hdlr:
  Handles ctrl-c signal in the threads
***************************************************************/
void sig_hdlr ()
{
   _exit(1);
}

/************************************************************************
                        DEF: matchends
  called by c2am.c from Ends->Ends button on menu
**************************************************************************/
void matchends()
{
  int i1, i2, k, j, c;
  CLONE *clp;
  struct match_info *tmp_match_list;
  int num_thread;
  int fd[500][2];
  int pid_array[500];
  int i=0,status, pid=-1;
  clock_t start=0,end=0, start2=0, end2=0;
  struct tms tmsstart,tmssend, tmsstart2, tmssend2;
  unsigned long num_entries;
  unsigned long actual_num_entries;
  unsigned int max_merge_status;
  unsigned int check_merge_status;

  sprintf(ZBuf,"\n>> Ends->Ends: Cutoff %.1e Auto %d Match %d  FromEnd %d\n",
          Pz.cutFlt, Pz.autoaddFlag, Pz.mergeMatches, Pz.fromendInt); Output;

  if (Pz.autoaddFlag) { /* cari 6aug04 */
     for(k=1;k<=max_contig;k++){
        if(contigs[k].count==0) continue;
        for(c=contigs[k].start; c!=-1;c=clp->next){
           clp = arrp(acedata, c, CLONE);
           clp->oldctg=clp->ctg;
           clp->ibc=0;
        }
     }
  }
  precompute_sulston();
  fix_cb_offsets();

  scanClam2();
  clearMsg();

  tctg = (struct tmpctg *) calloc(sizeof(struct tmpctg), max_contig+1);
  NOMEM2(tctg, "tctg");
  init_tctg();

  /* we will make a match list array that should be big enough */
  match_list = (struct match_info *) 
          calloc(sizeof(struct match_info), max_contig*MERGE_MAX);
  NOMEM2(match_list, "match_list");
  next_match_entry = match_list; 
  last_match_entry = &match_list[max_contig*MERGE_MAX - 1];
  tmp_match_list = match_list;
  num_thread=Num_threads;

  start=times(&tmsstart);
 
/** Parallel code **/
  if(num_thread > 0)
  {
   /* Calculate the number of comparisions to be performed per thread */
     update_end_to_end_struc(&num_thread);

     i=0;
     while(i<num_thread)
     {
        if(pipe(fd[i])<0)
        {
             printf("Pipe error \n");
        }
        if((pid=fork())>0)
        {
             close(fd[i][1]);
             pid_array[i]=pid;
             i++;
        }
        else if(pid<0)
        {
            printf("Fork :Error");
        }
        else
             break;
     }

   /* In the child */
     if(pid==0)
     {
        unsigned int thread_num=i,contig_num=0;
        unsigned int start_contig=end_to_end[thread_num].start_contig;
        unsigned int end_contig=end_to_end[thread_num].end_contig;
        struct match_info * match_list_iterator = match_list;
        unsigned int i2=0;
        num_entries=0;

      /* Need to register the signal handler for SIGINT (ctrl-C) 
             so that GTK does not crash*/
          signal(SIGINT,sig_hdlr);
      
          ZBufPrint = (char *)(malloc (sizeof(char) * ZBufSize ));
          ZBufPtr = ZBufPrint;
     
          for(contig_num=start_contig;contig_num<=end_contig; contig_num++)
          {
              for (i2 = contig_num + 1; i2 <= max_contig; i2++) 
              {
                  if (tctg[i2].n == 0) continue;
                  check_ctg_matches(contig_num,i2);
              }
              if (ZBufPtr != ZBufPrint)
              {
                 ZBufPtr=ZBufPrint;
                 if (Pz.logFlag) 
                 {
                    fprintf(Zlogfp,"%s",ZBufPrint); 
                    fflush(Zlogfp);
                 } 
                 if (!Zbatch_flag && Pz.stdFlag)
                 {
                    printf("%s",ZBufPrint); 
                    fflush(stdout);
                 }
                 memset(ZBufPrint,0,strlen(ZBufPrint));
              }
          }
               //Calculate number of entries in the match list
          while(match_list_iterator != next_match_entry)
          {
              num_entries ++;    
              match_list_iterator ++ ;
          } 
          close(fd[thread_num][0]);
          write(fd[thread_num][1],&num_entries,sizeof (unsigned long));
          write(fd[thread_num][1],match_list,
                    sizeof(struct match_info) * (num_entries));
          close(fd[thread_num][1]);
           _exit(0);
      }
      /* Parent Code */
      //wait for children
      //read from the pipe
      //update the common structure.

      i=0;
      num_entries=0;
      actual_num_entries = 0;
      max_merge_status = 0;
      check_merge_status = 0;
      while(i<num_thread)
      {
           num_entries=0;
           read(fd[i][0],&num_entries,sizeof(unsigned long));

      /* Check if it has exceeded max allowable merges */

           if (tmp_match_list + num_entries - 1> 
                        last_match_entry && max_merge_status == 0)
           {
              actual_num_entries = num_entries;

              //tmp_match_list points to the entry where 
              //the addition of new entry will start from
              // i.e if tmp_match_entry =1, last_match_entry = 3 , 
              // then (3-1) + 1 = 3 entries can be added
              num_entries = last_match_entry - tmp_match_list + 1;
              check_merge_status =1;
           }

           if(max_merge_status != 1)
           {
              read(fd[i][0],tmp_match_list,
                               sizeof(struct match_info)*(num_entries));
              if(actual_num_entries > 0)
              {
                 //Read the remaining data into a tmp buffer 
                 //It needs to be discarded
                 char * tmp_buffer;
                 tmp_buffer = (char *)(malloc (sizeof(char) * 
                                       sizeof (struct match_info) *
                      (actual_num_entries - num_entries)));      
                 read(fd[i][0],tmp_buffer,sizeof(struct match_info)*
                                (actual_num_entries - num_entries));
                 free(tmp_buffer);
                 tmp_buffer = NULL;
              }
      
              tmp_match_list = tmp_match_list + num_entries;

              // tmp_match_list can advnace one beyond last match entry
              if (tmp_match_list > last_match_entry)
              {
                 tmp_match_list = last_match_entry;
              }
              close(fd[i][0]);
           }
           else  /* Just read it so that the child exits */
           {
              char *tmp_buffer;   
              tmp_buffer = (char *)(malloc (sizeof(char) * 
                     sizeof (struct match_info) * (num_entries)));      
              read(fd[i][0],tmp_buffer,sizeof(struct match_info)*(num_entries));
              perror ("read");
              close(fd[i][0]);
              free(tmp_buffer);
           }

      /* wait for the child so that it doesn't show up as <defunct> */

           waitpid(pid_array[i],&status,0);

           if (WEXITSTATUS(status) == 1)
           {
           /* Child has exited because of ctrl-C being pressed */
           /* So kill all children and exit           */
              int j=0;
              for( j=0; j<num_thread; j ++ )
              {
                   kill (SIGTERM, pid_array[j]);
              }
              printf("Exit::User initiated abort") ;
              exit(0);
           }   

           if (check_merge_status && max_merge_status == 0)
           {
              sprintf(ZBuf,"Maximum number of matches exceeded\n");
              Output;
              max_merge_status = 1;
           }
           i++;
      }
      next_match_entry = tmp_match_list ;
      tmp_match_list = match_list;

      while(tmp_match_list != next_match_entry)
      {
          tctg[tmp_match_list->ctg1].num_merges ++;
          tctg[tmp_match_list->ctg2].num_merges ++;
          tmp_match_list ++;   
      }   
  }
  else /* Serial code - Number of threads is 1 */
  {
      ZBufPrint = (char *)(malloc (sizeof(char) * ZBufSize ));
      ZBufPtr = ZBufPrint;
      for (i1=1; i1 < max_contig; i1++)
      {
          if (graphInterruptCalled())
          {
              printf("User Interrupt - pre-mature termination of Ends->Ends\n");
              goto CLEANUP;
          }
          if (tctg[i1].n == 0) continue;

          for (i2 = i1 + 1; i2 <= max_contig; i2++)
          {
               if (tctg[i2].n == 0) continue;
               check_ctg_matches(i1,i2);
           }
           if (ZBufPtr != ZBufPrint)
           {
              ZBufPtr=ZBufPrint;
              if (Pz.logFlag) 
              {
                 fprintf(Zlogfp,"%s",ZBufPrint); 
                 fflush(Zlogfp);
              } 
              if (!Zbatch_flag && Pz.stdFlag)
              {
                 printf("%s",ZBufPrint); 
                 fflush(stdout);
              }
              memset(ZBufPrint,0,strlen(ZBufPrint));
          }
      }
  }
  if(ZBufPrint)
  {
      free(ZBufPrint);
      ZBufPrint = NULL;
      ZBufPtr = NULL;
  }
  end=times(&tmssend);
   
  for (k = 1; k <= max_contig; k++) 
  {
       contigs[k].projmsg[0] = 0;
  }


  if (Pz.autoaddFlag == 1)
  {
      /* set the oldctg parameter for each clone so viewer 
                can know where it came from */
      for (k = max_contig; k > 0; k--)  
      {
         for (j = contigs[k].start; j != -1; j = clp->next) 
         {
            clp = arrp(acedata, j, CLONE);
            clp->oldctg = clp->ctg;
         }
      }
      start2=times(&tmsstart2);
      qstuff = 1;
      do_automerges();
      qstuff = 0;
      end2=times(&tmssend2);
      fix_cb_offsets();
  }
  else
  {
      add_match_remarks(0);  
  }
  printf("\n");
  pr_end_timer("Calculate merges:", end-start,&tmsstart,&tmssend); 
  if (Pz.autoaddFlag) 
       pr_end_timer("Perform   merges:", end2-start2,&tmsstart2,&tmssend2); 

CLEANUP:
   free(tctg);
   free(match_list);
   Zproj_results(0); /* project window */
}
